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Abstract 


Convolutions on the sphere with corresponding cowolu- 
tion theorems are developed for one- and two-dimonsior?al 
functions. Some of these results are used in a study of 
isotropic smoothing operators or filters. Well known filters 
in Fourier spectral analysis, such as the rectangular, Gaussian, 
and Hanning filters, are adapted for data on a sphere. The 
low-pass filter most often used on gravity data is the rec- 
tangular (or Pellinen) filter. However, its spectrum has 
relatively large sidelobes; and therefore, this filter passes 
a considerable part of the upper end of the gravity spectrum. 

The spherical adaptations of the Gaussian and Hanning filters 
are more efficient in suppressing the high-frequency compon- 
ents of the gravity field since their frequency response 
functions are strongly tapered at the high frequencies with 
no, or small, sidelobes. Formulas are given for practical 
implementation of these "new" filters, including a demon- 
stration that the large negative sidelobe of the Pellinen 
response can cause 180 ® shifts in the smoothed function. 
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1, Introduction 


To the statistician the word ’’mean" denotes the value 
to be expected when sampling a given population of values. 
Statistically, the expected value is nothing more than one 
of several constant parameters that describe the population. 
When the geodesist speaks of a mean gravity anomaly (or mean 
geold undulation) he does not consider the entire terrestrial 
population of gravity anomalies, and yet it serves as a consit- 
uent descriptor of the earth's gravity field. The mean gravity 
anomaly is, in the most general terms, defined as the (possibly 
weighted) average of a subpopulation of anomalies distributed 
over a particular region of the earth's surface, for example, 
over a block delimited by pairs of latitude and longitude lines 
The total number of regions of a given size which together 
form the earth's surface is finite, but there exists an infin- 
ite number of ways to partition the surface into regions 
of one size (for example, by simply changing the location 
of the zero meridian). Consequently, the set of corresponding 
mean anomalies (the "moving average") forms a new infinite 
population which reflects the characteristics of the total 
gravity field to some degree of detail. It describes a field 
that* more or less, is a generalization of the actual field, 
representing the dominant or essential features and suppressing 
unnecessary or unwanted details. 

The following mathematical treatment sets the stage 
for the study of the different weighting schemes that can 
be used to define the mean gravity anomaly, thus making the 
above loose statements more rigorous. This can be accomplished 
effectively only by representing the gravity anomaly in terms 
of its spectrum, which is the set of coefficients in its 
representation by a series of spherical harmonic functions. 

This definition of the spectrum necessitates the approximation 
of the earth's surface by a sphere only if the spectra of 
the terrestrial anomaly and the anomaly on a sphere external 
to the earth (e.g. for satellite applications) are to be 
consistent. Otherwise, a spectrum is definable for any surface 
approximating the earth, that can be mapped onto the unit 
sphere viing a one-to-one correspondence. The gravity anomaly 
enters the discussion only as an example, since any of the 
geodetic quantities, indeed any function that is expandable 
in spherical harmonics, would serve equally well. 


2. Convolutions on the Sphere 


The first part of this paper shows how several concepts 
of spectral theory common in electrical engineering and com- 
munication theory can be applied and understood in physical 
geodesy (see also Robertson, 1978; Kaula, 1959, 1967). The 
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subsequent formulas by themselves are not new in geodesy » 
however they are viewed from the standpoint of spectral theory 
which acts to unify several diverse areas of physical geodesy. 

Since geodetic data pertaining to the gravity field 
are obtained either on the earth's surface (approximately 
spherical) or by earth-orbltlng satellites (approximately 
circular orbits), It becomes computationally expedient to 
consider functions in terms of spherical coordinates. Note, 
however, that for certain local studies the plane may serve 
as a sufficient approximation of the earth's surface. In 
which case two-or three-dimensional Cartesian coordinates 
are more appropriate (Moritz, 1966; Rayner, 1971; Breakwell, 
1979; see also Jordan (1978) who develops an Interesting 
synthesis of the global and local situations). The Cartesian 
coordinates are particularly' attractive due to the ease with 
which the Fourier transforms can be computed. 

Let F(8,X) be a function defined on the unit sphere 
(the scale of the sphere is immaterial). 6 , X are the 
usual spherical coordinates, being respectively the polar 
angle (colatitude) and the longitude. In analogy to the 
familiar Fourier transform, applicable to functions defined 
in rectangular coordinates, we define the two-dimensional 
'’Legendre transform" as 


L2[F] 


1 


// F(9,X) 


O 


^„„(e,x) 

nm 


da 


nm 


( 1 ) 


where a represents the unit sphere (0<X<2 tt, 0^9<tt), 
da * slnO dX dX , and where, for n>0 , 


Y (e,X) ■ P I 1 (cose) 
nm * ' n|m|' 


cosmX , 
sin|m| X, 


m > 0 
m < 0 


( 2 ) 


These are the (surface) spherical harmonic functions, which 
satisfy the following orthogonality property: 


1 


a 



if 

if 


n * 1 and m * k 
n a or m#k(3) 


In the mathematical literature the ?nm are sometimes defined 
using the exponential of a complex angle; however, the use of 
the sinusoidal functions as in (2) is more common in physical 
geodesy. The Pnm are the fully normalized associated Legendre 
functions: 
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(4) 


P„m<cos 0 ) P„„(cose) 

where er, *1 » em"i » w 0 > and where Pnm(cos 6 ) * sin *”0 • 

Hin ' 

*3TcoiF5^ Pu(cos0), and the P^ are the well known Legendre 

polynomials. The set of coefficients {fnm^ » definition, 
constitutes the (spherical Legendre) spectrum of F . As 
noted in the introduction, a different spectrum can be defined 
by replacing (not a formal change of variable) the geocentric 
colatitude by the complement of either the geodetic latitude 
or the reduced latitude. Since any function defined on the 
sphere is intrinsically periodic in both variables 9 and X , 
the spectrum is a discrete (but generally infinite) set of 
numbers (see Papoulis, 1976, p,70). 

The Inverse Legendre transform is then defined as 

» L - F(6.^) (5) 

n*Q m^-n 

The last equality follows only if F is continuous, but 
the series converges under less stringent conditions (Hobson, 
1965, p,342). Equations (1) and (5) express a duality between 
a function and its spectrum. Given the function, its spectrum 
is unique; and conversely, a given spectrum determines the 
function uniquely as long as the series (5) converges. Def- 
initions (1) and (5) differ from those of Robertson (1978), 
but the symmetry of the above transforms is too strong to 
resist. Unfortunately, this requires some sacrifice in symmetry 
for the one-dimensional Legendre transfo.rm and its inverse 
(m=0 in (1) and (5)); 

Li[F] = F( 0 ) p^(cos0) sin0 dO * f^ ( 6 ) 

go 

[f ] = I v'5n+I f P„(oos 6 ) = P(e) (7) 

n n n 

The coefficients (fn^ constitute the spectrum of F when 
the basis functions are the zero-order fully normalized as- 
sociated Legendre functions. Since a one-dimensional function 
defined on a circle is more conveniently transformed using 
the Fourier transform, the one-dimensional Legendre transform, 
per se, will not be considered further. Any function depending 
only on 0 is to be regarded here as a function defined 
on the sphere and Independent of X ; thus if 3F/3X = 0 , 
then F(0) = F(0,X) . Its spectrum, as given by (1), is 
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The spectrum of s function defined on a sphere depends 
on the orientation of the coordinate system. Consider the 
rotation of the coordinate system by the Euler angles a , 

13 p Y (see Fig. 1). Diverting momentarily to Cartesian 
coordinates x,y,z, the spherical harmonics are homogeneous 
polynomials (Kellogg* 1953f p.l39>* which under linear transfor- 
mations* such as rotations* transform into homogeneous poly- 
nomials of the same degree. The maximal set of independent 
spherical harmonics of degreo n comprises the 2n-fl harmonics 
of degree n and orders k * k«-n * ... * n . Hence the 
transformed spherical harmonic ^nmC'j*!^) the new coordinate 
system can be expressed as a linear combination of the harmonics 
^ nk^ ® I ^ ) » k*— n * ... * ns 


where the Cnit^^ are coefficients depending on the Euler 
angles (see also Muller* 1966; and Kaula* 1959). Cushing 
(1975* p.596) gives vise transformation coefficients explic- 
itly* but for spherical harmonics defined with the exponential 
of ccmiplex multiples of the longitude. Adapting this result 
to the definition of the spherical harmonics as given by 
(2)* it is only a matter of careful manipulation to derive 


•^^cm 


(- 1 )*" 

®nkm 



C? 4 

km 

(-1)*" 

c? 


«m 


- 

-km 

S + s" 

km 

km 

„n 


+ 

-km 

km 




■km^ 


k > 0 , 

k * 0 * 

k < 0 * 

k > 0 , 
k = 0 , 
k < 0 * 


m > 0 


m > 0 


m > 0 

( 9 ) 

m < 0 
m < 0 
m < 0 



where 


+ W) (0) 
®km * <*ton 


and 


dJm<6) - (-1)“""' 2™ '^{otTTThS^ oos8)*^a-oosS)T^ 


k-m 


n+m 


2 "*^ )l-0 


1 . <n;!li> <-l)"*'"-^l-oos0)'*"'-^l + oosB/ 


The sum in (11) including the. antecedent power of 2 is the 
Jacobi polynomial (cos|3) . Note that both 

integers k and m ”can assume positive as well as negative 
values and they must be taken literally in the context of 
each of the foregoing expressions (e.g. if k, m<0, then 
-k , -m are positive, and the argument of sln(kot+mY) is 
negative, etc.). Finally, we note that the binomial coef- 
ficients in the ur^pression for the iTacobi polynomial are 
defined such that (§) » 0 , if q>p , so that the summation 
extends only over those Indices A for which 


max(0, m+k) 1 A < min(n+k, n+m) 


( 12 ) 


for any k, m. If Y"0 » then 

<-!)“ -i- cosko I<-l>*‘dj| (3) +dV(l3)], k>0 , m k 0 


(-I)" 

♦'^m 


, k =* 0 , m > 0 


(-l)"* sin|k|a[(-l)*'d!]^^ffi) +dj||^(3)], k<0 , m > 0 


'nkm 




sinka[(-l)‘‘ d^<3)-d"^<3>) , k>0 , m<0 


0 


, k = 0 , m < 0 


coskot[-(-l)^d"^^(6) +d”^(3)] , k <0 , m <0 


( 13 ) 
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And ifoi-X, 0-5, T-0, m-0 (seo Figure 2), then we 

obtain the familiar addition theorem lor Legendre polynomials { 

P^(coe>l ) - 2n+i 


where cos’ll - cosO cos5 sinO sin5 cos(X-X). 

An important concept in signal and time series analysis, 
arising also frequently in geodesy, is the concept of con- 
volution. For example, the Stokes function convolved with 
the gravity anomaly results in the geoid undulation (see 
Robertson (1978) for other examples). The convolution of 
two one-dimensional functions F , G defined on the real 
line is formulated as 


H(x) - (G * F) (x) « / G(x-x) F(x) dx (15) 

-cc 


The convolution theorem in Fourier analysis is well known 
(Bath, 1974, p.79); it states that (aside frori a constant 
factor) the (Fourier) spectrum of the convolution H is 
the product of the (Fourier) spectra of F and G . The 
process of convolution in the space domain is thus transformed 
into the computationally simpler process of multiplication 
in the frequency domain. A similar result holds for convo- 
lutions on a sphere. 

Churchill and Dolph (1954) defined the convolution 
of two one-dimensional functions whose spectra are given 
by the Legendre transform. Since functions depending only 
on 6 here are to be regarded as two-dimensional functions 
on the unit sphere, but independent of X , the following 
definition of convolution differs from theirs; 


H 2Tr ir 

H(0) - (G*F) (9) » ^ f / G(t|;) F(9) sin0 d0 dv (16) 

tT] '0 0 

where cost|< - cos0 cos0 + sln0 sin0 cosv f. The variable v 
is part of the definition of the convolution and is necessary 
to derive a corresponding convolution theorem. Let , 

{g„} be the (Legendre) spectra of F and G , respectively. 
Then with the addition theorem (14) we have 


00 


G(V<) 


= I 


n 


g 


JL 


n-0 tOTTT m~-n 


l ?nm^e,v+J) 


nm 


(17) 
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for arbitrary % , In this and aXX foXXowing derivations 
the uniform convergence of the inverse transforms is presumed, 
so that the summation and integration may be freely inter- 
changed. Therefore 


00 n 

11(5) - I -In— I y (5.X) 
«i0 /5n+I m— n 


Sir Tr 

^ V-0 0-0 


(18) 


Substituting (2), the integral on the right side is zero 
unless m«0 ; hence with (6) 


m) 


f IQ-Ja 

n*0 /wn 


P ^ (cos5) 
no 


(19) 


The (Legendre) spectrum of H is therefore — f„ g„ , 

/2n+T 

not exactly the product of the spectra of P and G as 
in the case of the Fourier transforms. 

The convolution of a longitude-independent function 
G with a function F of both variables may be defined as 

- 2 tt it 

H(5,X) b (G * P) <5,X) » ^ /j G(tlj) F(6,X) sinO d0 dS 

( 20 ) 

Where in this case 


cosijj “ COS0 COS0 + sin0 sin5 cos(X-X) 


( 21 ) 


The corresponding convolution theorem is established by sub 
stituting (17) with = X - A into (20): 


H(5,X) 


oa 


I 

n=0 


g 


h— ? ^ (^,X) 

»^n+r m=-n 


1 


// F(0,X) da 


00 


n 

= l 1 

n-0 m*=-n 


gn 


•nm V 


/2n'+I 




(22) 


where (1) was used. The spectrum of H is therefore 
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(23) 


h 


nm 


v€iHT 


a result differing considerably from the one-dlmenslonal 
Fourier analogue. The analogy is almost completely lost 
because the spectrum of G is actually the set (gnm^ » 
which consists only of zonal coefficients, i.e. g«o = gn 
g„^ * 0 , if m # 0 . 

In the area of functional analysis, the convolution 
integral (20) is viewed as an operator, 

*’ = •A U 

with an associated kernel G . Thus rF = H . If the result 
of operating on a function is a scaled version of the function 
itself, then it is known as an eigenfunction; the scale 
factor is called the eigenvalue: 


rr = fF 


(25) 


Recognizing the set 
if i =yn or j m) 
function ?pm » 
(using (23)) 


{yij=l if i = n and j =m ; 
as the spectrum of the spherical 
spectrum of the convolution is 


harmonic 


1 _ 

/H+I 


Si 



/2H5T ” 

0 






(26) is the spectrum of the function 
ry = ^ — g ? 


if 

i != n 

and 

if 

i n 

or 

1 

/2inH 

i "" 

^nm 


3 = m 
j # m 

(26) 
so that 


(27) 


Therefore, from this slightly different perspective, the 
eigenfunctions and corresponding eigenvalues of the operator 
(24) are respectively the spherical harmonic functions 

coefficients g^/ /2n+T cf. Melssl, 1971). 

Finally, a similar definition holds for the convolution 
of two functions defined on the sphere, each depending on 
two variables; 


R(e,X) = (G*F) (e,X) = /J G(i/J,C) F(e,X) da ( 28 ) 

^ a 
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whore i}^i » C sphorieal coordinates » colatitude and 

longitude, in a system rotated by the angelos IT , T ; sin» 
Fig» 2. A comparatively simpl© relationship among the spectra 
ot two functions FCO ,X ) and G(0 ,X ) and their convolution 
(equations (28)) does not exist. I;C we substitute the trane- 
formation (8) with the Euler angles y*0 

into the spoctral representation of G , then (28) becomes 


«» n 


I I nOA) ^ 

a n«0 m*-n 


«* n 

n*0 m*~n 


K*-*n 


« n n 

nL) mw n k« n uT jS*) 

n*o m**<^n k"~n 


(29) 


The following is a formal derivation of the equation 
for the spectrum of H , oxcliidos a practical method of com- 
puting it, and can be omitted without loss of continuation. 
The final result is given by equations (37) through (tl), 
Considering (13), we may write 




X ) 



t coskX 
sin I k| X 


k >0 and m^ 0?, or, k< 0 and m< 0 

k <0 and m> 0; or, k> 0 and m< 0 

(30) 


whore the definition of »s{!m(S') can be inferred from (13). We 
note horo only that, according as ktm is either even or „ 
odd, n-th degree polynomial in cosff or sinO 

times a similar polynomial of degree n-1. This is verified 
by noting that the expression for d)',„(0') (equation (11), 
tt*0) contains only factors of the form (1 - cosO‘)*'"b • . 

(ItcosO)J^ (mtk even), or sinil (1 - cos^)h“P~a (i + cosli)^”® 

(m t k odd), whore p«A— J(ktm) li 0 (always). The ^nkm^^*^^ 
are therefore analytic functions and can be expanded as series 
of spherical harmonies; 

« i 




fio d»*-i 


(5,X) 


(31) 


By (1) and (5) 


i jnm 


II 

0 




/■coskX 
vsln (kl X 


) aSii) *“ 


> 10 * 


( 32 ) 


Invoking the orthogonality oi the sinusoids, we have 


C 4 4 * 

^ijnm 


0 J ,1 ^ k(m > 0) , d # -kCin < 0) 

(33) 

7T 

(cos5)slnQ‘d0 ; d-k(m>0), d«-k(in<0) 

(33) 

Hence, because i cannot be less than |kl , (31) becomes 


^nW0»^> * ^J|j^|^ik*nm 


(34) 


where k*»(sign of m)» k . Substituting this into (29) yields 

w n n “> 

5;„ ^ X , L ,«'nm ^nk ^lk»'°-^> 

nwQ m«-n k«=-n i« k 

(35) 

M n 

Now it is easily recognized that, symbolically, 1 X * 

CO 00 n*so !n«“n 

J y ; therefore, by first transposing summation signs, 
m=!-<» n“[m| 
we have 


oonnoo oon<^*” 

III I ^ I I II 

n«0 m=^-n k“~n i*[k| n«0 k»~n i*^| ms~n 

00 00 CO n 

1 I , I , I 

lca -00 n*=|k| i=|k| m«~n 

oo i 00 n 

II i , i 

i»0 k®-i n«|k| m-~n 


(36) 


It can then be verified that 

« i “ n 

■ .1 a. j., J-. »*«» 


(37) 


where 


^iknm 


(®nm *nk 'iknm > "> i 0 

8,™ *n,-k '^iknn, 
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( 38 ) 


' 'fr - 


and 


?lknm “ k| 


Finally I 

H(0,X’) 


W i 

“ I I 

i«0 k»-i 




(0 ,D 


where 




n 

J „ ^iknm 


(39) 


(40) 


(41) 


The set {hik> is the spectrum of the convolution (28). 

The evaluation of hji^j^ , given the spect.m of F and G , 
involves an infinite summation, as well as the computation 
of the integrals (39) . 

Clearly, from the above discussions, the powerful convolu- 
tion theorem of Fourier analysis cannot be adapted blindly 
to two-dimensional spherical functions. Although the spectrum 
of the spherical convolution also involves the products of 
the spectra of the functions being convolved, the relationships 
are generally not as straight forward as in the Cartesian 
case. 


To conclude this section the Dirac delta function is 
defined on the unit sphere. In the rectilinear case, this 
function, denoted 6(x) , is defined to be zero everywhere 
on the real line except at a single poiat x such that for 
a continuous function F(x) 

00 

/ 6(x-x) F(x) dx - F(x) (42) 

—00 

The spherical equivalent of 6(x) is defined similarly, 
but here we assume that it is also isotropic, i.e. indepen- 
dent of the direction between the point of integration and 
point where it is nonzero. Denoting it by D(t|j) , we have 
(by definition) 

^ // D(\|i) F(0,X) da = F(0,X) (43) 


where i|j is the spherical distance between (0,X) and (0,X) 
(see (21)). The spectrum of D(^) is (see (1)) 


12 - 


( 44 ) 




0 



D(i|;) P^^(cosip) da 


m 0 
m « 0 


Hence, using (43), 


dn = Pjj^(cosO°) = /2n+i 


(45) 


The operation associated with the delta function is a con- 
volution of the type (20). In the realm of functional analysis 
in Hilbert space, D(i|;) is an example of a reproducing kernel 
(see Krarup, 1969, p.43). 


3. The Mean Gravity Anomaly 


In this section some of the mathematical tools developed 
above are implemented to study the smoothing of the gravity 
field. The (point) gravity anomaly, with the standard spherical 
approximation already applied, has the following series rep- 
resentation (Heiskanen and Morit' , 1967, pp. 89, 108): 


00 

Ag(r,0,X) = y \ (n-1) (-^)”^^ j (Cj^ijjCosmX + S^i^sinmX )P^jjj(cos0 ) 


where r is the distance from the css>nter of the earth; 

Y = kM/R^ is an average value of gravity; kM is the product 
of the gravitational constant and the earth's mass; R is the 
radius of the sphere that approximates the earth's surface; 
and Cnm » Snm constant, dimensionless coefficients. 

Equation (46) may be rewritten in the more compact form as 


n 


Ag(r,9,X) = I Ann, 

n=o m=-n 


where the coefficients. 


Y(n - 1) Cnnd 
Y(n-l) 


ni > 0 
m < 0 


(47) 


( 48 ) 
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constitute the spectrum of the gravity anomaly on the sphere 
of radius R . It is worth noting that the spectrum on a 
sphere of radius r«Ri >R is given by the set of coefficients 

{(R/Ri)"+2 A„„) . 

From the introductory remarks, the mean gravity anomaly 
is derived from the point anomaly by subjecting the latter 
to an averaging process. Mathematically, we formulate this 
by applying an operator (only isotropic operators will be 
considered), such as, 

^ // B(t|))(- ) da (49) 


to the gravity anomalies within the region Aa on the earth's 
surface. Because the kernel of this operator, B(i|^) , is 
supposed to depend only on \p (i.e. the spherical distance 
between the point of computation and the point of integration), 
the region is necessarily a spherical cap centered at the 
point of computation. By defining B(i|i)s=0 outside the 
cap we obtain the following more general formulation of the 
weighted average; 

Ag(0.X) ^ B(ip) Ag(0,X) da (50) 


The operator thereby becomes a member of the class of oper- 
ators (24), and the operation itself is a convolution of 
the type (20). For convenience, we may write the kernel 
as a normalized weighting function; 

B(l|0 = -r, (61) 

Ji // W(I)I) da 

a 

where |w(i|j)| for 0 < ij; < ir . Now let {bjj} be its spec- 

trum, i.e. 


00 

- I /2n+T bj^ Pjj(cosif)) (52) 

n=0 


where, according to (6) and (7), 


b 


n 


/2nTT 

”"2 


/. 


B(>|0 P^(cosiji) 


simj; dip 


(53) 


The convolution theorem, enunciated by equations (22) and 
(23), then directly provides the spectrum of the mean anomaly 


Ag ; 

A 

nm 


( 54 ) 


are the eigenvalues of the averaging oper- 
ator; or in the jargon of spectral theory, the frequency 
response of the filter. In view of (54), | 100 0n| is the 

percentage of an (n)th-degree harmonic coefficient that is 
retained in the process of averaging. 

Averaging gravity anomalies over a spherical cap is 
not a typical operation in geodetic practice. Rather, mean 
anomalies are viewed as averages over spherical blocks (trape- 
zoids) delimited by convenient global or local curvilinear 
coordinate lines. The average over such blocks is characterized 
by operators that are nonisotropic (see Pollack, 1973), and 
the averaging process is then a convolution of the general 
type (28). The spectrum of such a mean is difficult to obtain 
in terms of the point anomaly spectrum (see the previous 
section). Gaposchkin (1980) derived and studied a series 
expansion of the mean over a trapezoidal block (however, 
it is not a spectral representation as defined by (5)). Pollack 
(1973) compared fche 3n corresponding to isotropic operators 
with the degre variances of the above mentioned nonisotropic 
operators. Recalling (41), it appears uncertain whether 
a study of these degree variances is indicative of how the 
spectrum of the anomaly is transformed in the averaging process. 
When computing mean anomalies from the point anomaly spectrum, 
the average over a block is frequently approximated by an 
average over a cap having the area of the block (Rapp, 1977). 
This approximation degrades with the elongation of the block 
(e.g. in the polar regions, where the meridional coordinate 
lines converge). Some further study is indicated here, but 
is beyond the scope of this paper. 


'n 


nm 


/5rvfl 


*= B A 

^n nm 


The constants 3n 


3.1 The Pellinen Mean 


The equally weighted average of the gravity anomaly 
over a spherical cap is defined as the integral of the anomaly 
over the cap divided by the cap's area. For this simple 
average (see also Pellinen, 1966), we have, in agreement 
with (50) and (51), 


1 


Wp(l|^) 


0 


0 < t|; i i|;o 

Vo ^ 
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(55) 








i|>o being the radius (generating angle) of the cap. Substi 
tuting the resulting kernel, 




4TT 


2rr sinij; d4» 


(56) 


into (53) yields 


.^0 


/ On4.1 -TV 

b » / P^(cos'P) sin'll 

1- cos'll 0 “ 


Pn-l(cosi|/6 ) ■ Pn-fl (cosi|;o ) 

/2n+l (1 cosipo ) 


(57) 


for which Sjoberg (1979) has found the following recursion 
formula; 


3p^ = 


n 


^Pn ^ 2n-l 


/2n+I 


n-2 


n+r ^P„_1 “ n+r ®P„_2 , n i 2 


Sp( 


1 » 3pi = 4(1 + costjio ) 


(58) 


Fig. 3 depicts wp as a function of \p for ijro =0?564 , 
which corresponds to a cap whose area equals the area of 
a l®x 1° block at the equator. The coefficients 3p„ are 
shown in Fig. 4. Here (as in Fig. 5), for the sake"of clarity, 
the values of the smoothing factors 3n (which are discon- 
tinuous functions defined only for integer values of n ) 
are traced by smooth curves. 

The spectrum of the Pellinen mean has an infinity of 
components, but the high-degree coefficients are clearly 
diminished in magnitude and thus suppress the local structure 
of the anomaly field. The first zero in the Fourier spectrum 
of the analogous one-dimensional "rectangular" filter occurs 
at n = = 319 for t|;o= 0?546 (Bath, 1974, p.218). This 

formula is not applicable in the spherical case as 3pn^*l^o ” 
0?564) is first negative when n = 389 . The major drawback 
of this operator, viewed as a "low-pass" filter, is the presence 
of the large positive and negative "side lobes" in its spectrum. 
That is, the filter admits, or passes, a considerable part 
of the upper end of the spectrum, even changing the sign 
of some of the coefficients. This "reversal of polarity" 
could transform minima of the point function into (false) 
maxima of the mean function, and vice versa (Holloway, 1958); 
see section 4. 


3.2 The Gaussian Mean 


This and the following weighting function are adaptations 
of filters commonly found in electrical engineering (Harris, 
1978, gives a good summ'^ry). The Gaussian mean takes its 
name from the bell<-shaped normal (Gaussian) probability density 
function that its weighting function resembles for small v ; 


WQ(t|j) - “ costfJ) ^ a>0, 0<i|i<iT (59) 

« e” T , small ijj 


The dimensionless parameters "a" (identifiable with the inverse 
of the second moment 0^ of the Gaussian distribution) charac- 
terizes the smoothing process. Since 



the kernel 
Bq(^> 


Wq(i]j) da s ^ (1 - e-^^) 

becomes 

2ae-a(l- cos»|^) 


(60) 


(61) 


and the eigenvalues of the corresponding operator are 


'n 


= *’Gn = 
/2n+T 


ae 


-a(l-y) 


-1 


(1 - e 


P^(y) dy , y^ cost|) 


(62) 


for which the following recursion formula holds: 


^Gn+1 
^Go = ^ 


2n+l 


'^n '^n-1 

-2a 


ft = ^ Q 


a 


(63) 


(63) is readily verified by applying the relationship P^+^(y) = 
(2n+l) Pn^y ^ “ Pn-l^y ^ ^®2) and integrating by parts. 

The function (59) is a global weighting function to be convolved 
with gravity anomalies over the entire sphere. On the other 
hand, restricting the averaging process to a cap, 
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^ [degrees] 



(64) 


Wq(i|;) 


Wq(i|^) 

0 


0 < 

< TT , 


leads to eigenvalue© satisfying the following recm’sion formu- 
la (see Appendix A)! 


'Cn'M 


0 * ^ 


2n+l 

a 



fl 1 « . 

®Gii ■*■ ■*■ j^_g-aa-y,, 

_ X-yaO-“<i-y'> 1 . 

l.;-‘a-y.> - T - 


n > 0 

yo •* costl^o (65) 


""he eigenvalues 0Gn (equation (62)) are all positive. 
To pruve this it is enough to show that 


I„ * e^y P„(y) dy>0; n>0, a>0 (66) 

11 n 


Substituting the uniformly convergent series expansion of 
the exponential function, ^ 1 ^ integral 

above yields k»0 

" in ■^1 t <-!>“■'" C y*' ^<>'1 “S' > 

la 1^ 

where we note that by orthogonality 


Ii dy « 0 , k * 0,..., n-1 


(68) 


and that y^Pn(y) is an odd or even function according as 
k + n is odd or even. Now, Hobson (1965, p.40) derived 


C ^n‘5'’ “y “ 


k(k-lKk-2)...(k-nt2) 

( k+n+1 ) ( k+n-1 ) . . . ( k-n+3 ) 


(69) 


This is always positive for k > n ; hence all terms of the 
sum (67) are nonnegative, thus establishing the inequality 
(66 ) 0 The signs of the harmonic coefficients are consequently 
preserved in the process of smoothing. 

The same is not true for the weighting function (64) 
as demonstrated in Fig. 5. The four cases shown here represent 
several possibilities to choose the parameter "a". With 
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1{»0 *0?564 , each of the conditions wg(iPo)»0.5 , 

0.1 , and 3a » \Ijq (where a «l/a*) produces different fre- 
quency responses (cases II, III, IV). Also shown (case I) 
for comparison are the 0 g- corresponding to the global 
weighting function with the condition wg( 0?564) »« 0.5 . The 
weighting functions associated with these smoothing factors 
are shown in Fig. 6. 

Figures 5 and 6 suggest that as the cap-edge value of 
the weighting function Wet'll) approaches zero (it can never 
equal zero), the oscillations of the correspondinf^ frequency 
response decrease in magnitude. Furthermore, the smaller 
"a" is, the more the higher-degree harmonics are filtered 
from the anomaly. Designing the Gaussian filt^,r so as to 
have ^artain smoothing properties is therefore accomplished 
by properly choosing values of “a" and 'po . The selection 
of '‘a*’ controls the essential bandwidth of the filter, i.e. 
what frequencies should be passed; while the subsequent 
choice of t{>o controls the magnitude of the oscillations 
of the frequency response function. The frequency response 
of the Gaussian filter defined for data on the real line 
is given by Holloway (1958, p,359) as 



(where, if f denotes frequency, the relationship 2iTf»n 
was used). The similarity between $Gn (global filter) 
and is shown in Table 1 and can be invoked to provide 

an approximate value of "a” given a desired value of 3n 
for some n . For example, if the harmonic at degree n 
is to be suppressed to lOOfoX of its original value (less 
than lOOfo/i of each subsequent harmonic will be passed), 
then 


“ T^HTTro 


(7X) 


Table 1; Legendre vs. Fourier Gaussian Frequency Responses 


■ 

a = 128235 1 

— 


a - 13131 


n 

5g„ ■■ 

(eau.(63)) 

(eau.(70) 


n 

BG„ 

(equ. (63)) 

Sn 

( equ . '(TO),), 

10 

.9996 

.9996 


10 

.9958 

.9962 

50 

.9901 

,9903 


50 

.9075 

.9092 

200 

.8549 

.8556 


75 

.8049 

.8072 

400 

.5350 

.5359 


100 

.6807 

.6833 

600 

.2451 

. 2457 


125 

.5490 

.5516 

700 

.1476 

.1480 


150 

.4221 

.4245 

800 

.08220 

.08246 


200 

. 2164 

.2180 

900 

.04235 

.04250 


250 

.09168 

.09256 

1000 

.02018 

.02026 


350 

.009300 

.009424 
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0.3 Hanning Mean 


The weighting function 




|i(l + cos«J|-) , 
0 


0 1 1 'I'O 

\|»o < i|j < ir 


(72) 


resembles the Hanning filter whose desirable filtering charac- 
teristics (Holloway, 1958; Bath, 1974, p.l58) make it popular 
in the disciplines where it is applicable. (Sen also 
and Arabelos (1981) where the Hanning "window function" (72) 

Is used to reduce the integration ("truncation") error incurred 
when performing an integration of a covariance function over 
loss than tho required interval, 0 ^ i|) < ir . The spectrum 
of wn defined by them differs essentially from the definition 
adopted here, but the spectral properties are similar.) 

The eigenvalues of the corresponding Hanning smoothing 
operator are given by 


ir 


2 


_ 

'll. It2(l -COSl|lo)-^'l''o^ 


1 


n 


where 


Ho 


*! 3L 


(73) 


P^^(cosi|/) simj} di|j , b * (74) 

A recursion formula for an (bs^ integer) is derived in Appendix 
B; it is 


[(n-M)2-b2 ]a,j-[(n-2)2~b2 3«n^2* Cntl)-Dj^-(n-2)“l)j^^2 


where 




n > 1 


» sinb iji 0 s ini|j o Pj^ ( eosift q ) 


, nil 


(76) 


1 


Co [l-GOS(jjo cosbt|io -bsiinjjo sinbipo] 

Cl ® 27T^') C2-2cos2i{ic) Gosbif^o *“ bsin2i}^o sinbi|)c] (77) 


a 


2 * [3 -3cos3'j*o cosbi'o - bsin3'jJo sinb>j'o] 

•~ 24 *' 


i 


a 0 


The first zero of the Fourier spectrum of the Hanning 
filter occurs at n - 2 Tr/<|^o « 638 , for i|io - 0?564 . For the 
spherical adaptation, BUnC'l^o wq? 564) first changes sign 
at n * 695 . Thereafter, the side lobes are comparatively 
small. Figures 4 and 3 show respectively the frequency response 
and the corresponding weighting function w^ . 


3.4 The Upward Continuation (Poisson) Operator 


Functions harmonic in the space external to a sphere 
centered at the coordinate origin attenuate radially according 
to r"Cn+l) . rAg is such an harmonic function (see equation 
(46)), and the gravity anomaly field at satellite altitudes, 
for example, is much smoother than at the earth's surface. 
Therefore, for a fixed radius, r * Ex , the upwa^"d continuation 
operator acts like a atmoothing operator with eigenvalues 
(recall the comment after (48)) 


3 « , n > 0 (78) 

Un Ux 


For the disturbing potential the exponent is n + 1 , and for 
the geoid undulation it is n - 1 ; in all cases Buo I , 
so that, unlike the previous filters, this smoothing process 
does not leave the global average unaltered (unless Aoo =0) 
The corresponding smoothing kernel is 

00 

B^Ojj) = (2n+l) (^)‘^'^^ P^(cosi|i) , 0 < < IT 


s^(i - ) 

[l - 2SC0S1|; + S2]V2 


(79) 


where s = R/Rx and the familiar generating function for 
Legendre polynomials was applied. (79) is nothing but s 
times the kernel of Poisson's integral (Heiskanen and Moritz, 
1967, p.35). Because the smoothing is biased (the integral 
of 4ir Bu(i|)) over the unit sphere is not unity), a weighting 
function satisfying equation (51) cannot be defined. However, 
for purposes of comparison, let 


w^(ij;) - 
Solving for 



_ (1-s)^ 

b“(0) [I- 2s cos 4- sf % 

s , we have 

4 (l|j) COSljJ - 1 + (2wJ^ 3 - cos^ji) - W^'4 

(ij;) -1)”^ 


(SO) 

(t|f) Sin2i|,)^4 
(81) 
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Although the ontiro gravity anomaly field must, in theory, 
bo convolved with the kernel (79) to yield the smoothed field, 
a specification such as w^(^{'o)“.01 renders the operation 
practically a local one. With >|i * iliQ •* 0?564 , eciuation (81) 
gives s«0.99785>0813 (or if R « 6371000 m , then Rf 6384851 m.) 
w„ and 8un with this value of s are shown in Figures 3 
and 4. 

The upward continuation produces a smoothed field with 
a more substantial attenuation of the low-degree harmonies, 
but with a fairly weak taper in the high end of the spectrum. 
Therefore, it is essentially different from, and less useful 
for smoothing than, the more conventioiml (weighted) areal 
averages, 


3,5 The Ideal Filter 


Ideally, the smoothing process should suppress the higher- 
degree harmonics completely with perfect retention of the 
lower-degree harmonies. This type of smoothing is automatically 
implemented whenever a determination of the gravity field 
consists of finding, through satellite techniques, the coef- 
ficients of its series representation. Obviously, only a 
finite number of coefficients can be found and the corresponding 
truncated series represents a perfectly filtered field (neg- 
lecting aliasivig and measurement errors). The eigenvalues 
of the operator that produces such ideal smoothing arc 


1 , n » 0 , , . , , n 

0 , n > n (82) 


where n is the desired degree of truncation — the larger 
n is, the less smooth is the filtered field. If the operator 
is global, the kernel is given by (see (52) and (54)) 

n 

B.(ijj) * I (2n+l) P„(cosi|0 (83) 

^ n=0 " 


(If n ^ 00 , Bj becomes the spherical equivalent of the 
Dirac delta function^ see (45).) 

It is not possible to define a kernel whose convolution 
with gravity anomalies in a cap with radius 0<t|Jo < it yields 
a perfectly filtered field. That is, a kernel with a finite 
spectrum is analytic on the sphere and cannot be gioro over 
any part of the sphere unless it is zero everywhere. 


3.6 Discrete Operators 


The averaging process in actual practice deals with 
discrete values of the gravity field — integrations such as 
(50) are performed numerically often using simple midpoint 
formulas: 


A"g(0 ,X) 


1 


M 

I B(ijjj) Ag(0 j ,Xj) ACj 

i 1 **1 


( 84 ) 


where 'j/j is the geocentric angle hetween points (S',X) and 
(0-j,Xj) , (Oj,Xj) is the center point of the area element 
Aoj > and M is the total number of Aoj into which the 
cap Oq has been partitioned. The above ‘^formula can be 
rewritten as an integral by utilizing the Dirac delta function 
D(i|^) (see (43)): 


M 


Ag(0,X)= {^)l B(tjJj) n D(Jj) Ag(0,X) da A0J 
= 


= (^)// ag<6,X) do 


(85) 


where cosij^ j - cos0 cos0j + sin0 sinOj cos(X-Xj) and iji , 5 

are the coordinates of the point (0,X) in the* spherical coor- 
dinate system whose pole is the point (^,X) (see Fig. 7). 


In Fig. 7, (Qj»Xj), M are fixed points and 

as the variable point (0,X) moves over the sphere, its position 
with respect to the point of computation (0,X) may be described 
by the coordinates » whence the arguments for the 

kernel B . Note that B is zero unless (0,X) coincides 
with one of the points (‘t>j>Xj). Thus the operator is noniso- 
tropic and the convolution is of the type (28). This implies 
that the frequency response of any of the filters discussed 
in the previous sections, in practice, is not given by the 
corresponding 3n • However, if the numerical integration 
is sufficiently accurate, the general characteristics of 
the response are retained, as seen in the next section. 
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pole of 
<0,X) system 


(0,X) 


Figure 7 ; The relationship of (0,X) to (0 j,X^) 
and equation (85). 


4 • Smoothing a Simulated Gravity Field 



In this section the formulas for the various averaging 
processes are put into practice, illustrating, by an example, 
the differences between the Pellinen, Gaussian, and Hanning 
averages in the space domain. For this purpose, a gravity 
anomaly field was generated as a series of spherical harmonic 
functions with coefficients; 




[(180,180) solution of Rapp (1978)], n=2,...,100 
"nm /■a“ > "=101,..., 360 (86) 


where the Unm are random numbers uniformly distributed on 
the interval [-0.5, 0.5] and d^ is their degree variance. 
Cji is the degree variance of tne gravity anomaly and is 
modeled (Rapp, 1979) by 
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- 3.405 ( .998006)"*^ + 140.03 ( .914232)"'^® 

tnigal]“ , n > 3 (87) 


The (X80,180) solution was derived from l^x 1° mean gravity 
anomalies, and the division hy 3pj^ , being the frequency 
response of the Pellinen averaging operator with ipo «0?564, 
effects a relatively smooth transition frop the actual (orig- 
inally mean) spectrum to the modeled (point) spectrum. The 
smoothing of the coefficients by a Gaussian filter (a® 14000, 
i|^o “tt) attenuates the upper spectrum to virtually zero at 
n « 360 ( 3 q 3 6 0 « 0.01) , resulting in a less irregular field 
for which the differences in the filters are better illustrated. 


The gravity anomaly function defined by (86) was evaluated 
on two profiles near the equator*; , 0?5 < A < 42S5 

(see Figures 10 and 11). This function was convolved with 
the Pellinen smoothing kernel (ijJo -1?692, corresponding to 
a cap having the area of a 3°x 3° block at the equator), the 
Gaussian kernel (ij^o =2?459, a = 4887.27), and the Hanning 
kernel (i|io = 2?459) . These latter choices for the cap radius 
and the parameter "a" yield frequency responses that pass 
generally the same band of frequencies as the Pellinen average; 
see Figure 8. The precise vfc,lue of the cap radius was estab- 
lished by requiring the cap to have the area of a finite 
number of 0?5 x 0?5 blocks arranged roughly in the shape 
of a disk, as in Figure 9. The parameter "a" then follows 
from the stipulated relationship 3a = i|;o (a^=l/a). 


The mean gravity anomaly function, computed using the 
smoothed coefficients 3nAnm' ^t’e plotted in Figures 10 
and 11 (Pellinen vs. Gaussian averages) and Figures 12 and 
13 (Pellinen vs. Hanning averages). The essential feature 
of these comparison is the 180° phase shift of the Pellinen 
average with respect to the point function at a wavelength 
of about 2° , as predicted by the frequency response 3p„ 

(in Figure 8, 3p„ has a maximum negative value at n = I72, 

corresponding to a wavelength of about 360/n = 2?l). The 
Gaussian and Hanning mean functions, on the other handj repro- 
duce more faithfully the peaks and valleys of the point function, 
as the corresponding frequency responses have either small 
or insignificant negative values. 


The averaging process, in the practical situation, is 
performed on discrete values of the point function. Consider 
the following discrete smoothing kernels 


, k = P, G, H 

♦ (where cJ) = 9O°-0) 


( 88 ) 
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Spherical caps of radii i]jq= 1?692 and i|Jo= 2°459 
defined as having the area of 36 and 76 0?5x0?5 
0?5x0?5 blocks, respectively. For the 
simulation, the continuous function is averaged 
over the caps, and the discrete function is 
averaged over the block values. 


cosijjj = sini^ sincjjj + cos^ cosc{)j cos(x-Xj) 
'^J+1 ’ " \l+l - \l 
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original (point) anomaly 



[degrees] 












whero 




(89) 


1 M 




and A 0 j is the area of the block containing a known value 
of the point function. is the discretized integral 

of the weighting function over the cap (see (51)); and although, 
for the Wi^(ij^) considered here, these Integrals are evaluable 
analytically, their discretization guarantees the unbiasedness 
of the averaging process 


Jn order to simulate the process of averaging, the coeff- 
flcients (86) were synthesized on a 0?5x0?5 grid in the 
equatorial region defined by -3?75 < < 3?75 , and 

0?25 < X s 42? 75. The Pellinen, Gaussian, and Hanning smoothing 
kernels were evaluated at ip j , j » 1, . . . , Mj^ , where 
is the number of 0?5x0?5 blocks in each of the caps (see 
Figure 9, Mp»36, Mg*M|j“ 76). Subsequently, the average 

gravity anomaly for each of the weighting functions (P,G,H) 
was computed on the profiles ^*+l?5, at 0?5 intervals, 
using equation (84). These "moving averages" are plotted 
against the original point function in Figures 14 and 15 
(Pellinen vs. Gaussian averages) and Figures 16 and 17 (Pellinen 
vs. Hanning averages). We note that, although the frequency 
responses of the (discrete) operators are not precisely those 
of Figure 8, the averaging characteristics illustrated in 
these figures are essentially identical to those described 
for the continuous (isotropic) operators. 

Normally, the equally weighted (Pellinen) average, in 
the example above, would be computed on a 3°x 3° grid. The 
"polarity reversals " would then not be directly apparent. 

Also, if the errors of the original point data were statis- 
tically uncorrelated, so would be the errors in the 3° Pellinen 
means. The propaged errors in the Gaussian and Hanning averages, 
however, will be correlated if the caps are not disjoint. 

Figure 18 shows the possible configurations for the above 
example if the averages are computed on a 3®x3®grid. Each 
shaded region represents the overlap between two "caps" which 
is responsible for the correlation. Assuming that all point 
data have the same standard error o with zero correlation, 
the correlation in the means is given by 

P^^(A,B) « (90) 

^ 0j^(A) 0j^(B) 
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where Lar is the number of points ji common to the two 
"caps", A and B (in Figure 18, » LAC“20); and 


q^CA.j) = ^ 40 j 


= "ifd'ftpAy.i 
4JrJ^k(Vo ) 


( 91 ) 


tjjA-j being the central angle between the center of the "cap" 
A and the point j ; and where 


, M 

aic(A) ® cJ^ J qj,(A,j) (92) 

j«l ^ 


If Aoj and constant, then 


Pj, (A,B) 


Lab 




With respect to Figure 18 
Pq(A,B) = 3.0 X 10"'' 
Pq(A,C) = 3.2 X 10“2 


and the above expajnple 

, pjj(A,B) = 2.3 X 10"" 

, Pjj(A,C) = 4.6 X 10“^ 


(93) 


(94) 


These are generally not considered significant correlations. 


5. Conclusion 


From its birthplace in electrical engineering, the appli- 
cation of the spectral theory has spread to many unrelated 
disciplines, becoming indispensable in the methodical analysis 
of large amounts of data. In geodesy, while there is the 
complication of having to work on the sphere, many of the 
basic concepts of spectral theory are directly applicable, 
as in the study of smoothing operators. The effects of the 
different methods to smooth the gravity anomaly field can 
be properly assessed only by inspecting the spectrum of the 
corresponding weighting function. It is shown that the mean 
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gravity anomaly acquires new and different connotations when 
various weighting schemes arc introduced. By an example, 
we saw that the excursions of the Pellinen frequency response 
into the negative can cause 180® phase shifts in the smoothed 
function (see Figures 10-17). Furthermore, from Figures 4, 5 
and 8 it is readily evident that both the Hanning and Gaussian 
smoothing operators (with a suitable choice of the parameter 
•'a*') do a better job of filtering the high-degree components 
from the total field than the Pellinen operator. Finally, 
it is noted that a truncated spherical harmonic series (e.g. 
nmax®^®®) cannot justifiably be characterized as a Pellinen, 
or other, mean function, unless the coefficients of the series 
are multiplied by the corresponding frequency response which, 
in addition, must have no significant components beyond degree 

”max • 
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Appendix A 


Let 
y « 

Vo " 


Substituting the relationship 


(2n+l) P^(y) « Pi^;i^(y)- Pi_i(y) 


where the primes denote differentiation with respect to 
into (Ail) and integrating by parts results in 


“ -2i^ - “(^n+1 ' ^^n-1 


Equation (65) follows upon realizing that, with y ~ cosi/f 
Vo = cosil^o , 




G„ “ ^ “ 7 . -a<l-y> P„(y)dy 

/2ri+i e y^dy yo 

yo 


ae 


-a 


1 _ 
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(Ail) 

(A. 2) 
y I 

)] 

(A. 3) 


(Ai4) 


Appendix B 


Let 

» f cosbtj^ P„(cos<jj) sin'/' d<jj (B.l) 

n o n 


Substituting (A. 2) into (B.l) and Integrating by parts yields 


“n “ "’ll ^^n-1 ■ ®n+l> 


where 

I’n “ [P„.l<oos*o) - P„^^(oos^.)l (B.3) 


and 


6n sinbi/^ Pj^(cost|j) sini|> d\|j 


(B.4) 


Now from (Hobson^ 1965, pp. 32-33), 

i(l-y“) P'(y) = -5^ (P„_i(y) - P„+i(y)> 


(B,5) 


Hence 

■^<^-1 - = T sln^ ^ P„(oos*)<H, (B.6) 


Again, integrating by parts, we obtain 


2hS «n-l - «n+l> = 

1 1 b 

■* TT / sinbij^ cosiji P^Ccosip) sinijj di|; + -- a„ ,(B.7) 

nxiii'Q n tin 

which, upon substituting (n+1) Pn+i(y) + ”^n-l^y^ ~ (2n+l)yPj^(y ) , 
as well as (B.2), gives after several manipulations 


(n+1)^ - b^ X. 
5nTT °n+l 


n^- b^ 
5n+T 





(B.8) 


where 

E„ = sinb^/o sin*, P„(oos*,) 
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(B.9) 


Equation (B.8) can be reformulated as 


«n-l - ^n+1 - * »“n> 

Combining this with (B.2), we have 

(n* -b*)a„ - n=D„ + - bE„ 

and 

(<n-2)^ -b2]a^^2 " (n-2)'^ " ^V2 


(B,X0) 


(B.Xl) 


(B,X2) 


The fi^’s are eliminated by subtracting (B.X2) from (B.XX) 
and substituting (B.XO), This results in 


[(n+l)^' -bMa„ “ C(«-2)' -b*]a „_2 + (n+1)* D„ - (n-2)“ D„_g + 

-b(Ej| -E jj_2. (B.13) 


thus proving equation (75). 


